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ABSTRACT 


The  implementation  of  QR-algorithm  on  ILLIAC  IV  is  described.   An 
ASK  subroutine  for  computing  all  eigenvalues  of  a  real  Hessenberg  matrix  of 
order  less  than  or  equal  to  6k   by  this  algorithm  is  attached.   The  QR-trans- 
formation  consists  of  the  decomposition  of  the  matrix  A^  into  the  product  of 
a  unitary  matrix  Q,  and  an  upper  triangular  matrix  R  ,    and  forming  A,     by 
post-multiplying  R  by  Q,  ,   where  A.  =  A  is  the  original  matrix.   All  eigen- 
values are  either  isolated  on  the  diagonal  or  are  eigenvalues  of  a  2  x  2 
diagonal  submatrix  as  k->  «>  . 
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INTRODUCTION 


An  ASK  program,  for  finding  all  eigenvalues  of  a  real  Hessenberg 
matrix  of  order  less  than  or  equal  to  6k   is  written  and  tested  on  the  B5500 
simulator.   The  QR- algorithm  for  real  Hessenberg  matrices  described  herein 
is  that  of  Martin  et  al  [2].   The  ILLIAC  IV  computer  time  required  for  per- 
forming one  iteration  (computation  of  A    from.  A,  ,  refer  to  k.l)   on  a  6k   x 
6k   matrix  is  approximately  10  millisecond. 

The  necessary  information  for  using  this  program,  is  given  in 
Section  2.   The  test  result  of  this  program  on  a  k   x  k   real  Hessenberg  matrix 
is  given  in  Section  3«   In  Section  k,    the  outline  of  the  QR-algorithm  is 
given  and  Section  5  is  devoted  to  the  actual  programming  technique  to  imple- 
ment this  algorithm  on  the  ILLIAC  IV  computer.   The  flow  chart  and  ASK  program 
are  attached  as  Appendices  1  and  2  respectively. 


2.   USAGE 


This  program  assumes  that  the  given  real  Hessenberg  matrix  is  stored 
in  the  core  memory  in  the  straight  storage  scheme  so  that  each  row  is  stored 
across  the  PE's,  starting  with  PEO.   The  real  and  imaginary  parts  of  the  eigen- 
values found  are  to  be  stored  in  the  two  rows  in  the  PE  memory  specified  by 
the  user.   The  original  matrix  is  destroyed  and  replaced  by  the  matrix  which 
results  from  the  QR  transformations.   The  content  of  ACAR2  and  ACAR3  are 
destroyed  since  the  ACAR3  is  used  for  linkage  between  the  subroutine  and  main 
program  and  ACAR2  is  used  for  passing  the  address  of  parameters  to  the 
subroutine. 

2.1  Calling  sequence 

Calling  sequence  for  this  subroutine  is: 

CALL  HOP  (N,  A,  WR,  WI,  IT) 
where  A  designates  the  first  row  of  the  matrix  and  is  declared  in  the  main 
program  as 

A:   DATA       aQ0,  aQ1  , ,    a^^,  (0.0)M, 

ai0>  ail  ' '   al,H-l,  (0'°)M' 


Yi,o'  Vi,i'  '  •  '  '    Vi,n-i'(0,0)M; 

where  M  =  6k   -  N. 
N  is  the  size  of  the  matrix  declared  as 

N:  EQU  6k; 
or 

DEFINE  N  =  6k  ##; 
or  given  as  an  integer,  i.e.,  6k. 

WR  and  WI  are  the  rows  in  the  PE  memory  where  resulting  real  and  imaginary 
parts  of  the  eigenvalues  respectively  are  to  be  stored.   They  are  declared  as 

WR:   BLK  1;   and  WI:   BLK  1; 
IT  is  the  row  vector  in  memory  to  which  the  number  of  iterations  required  for 
finding  each  eigenvalue  is  to  be  placed.   This  is  declared  as 

IT:   BLK  1; 


3 
If  zero  is  placed  in  place  of  IT,  it  is  to  be  considered  that  a  user  does  not 
want  to  know  the  number  of  iterations  required.   The  CALL  macro  should  be 
defined  in  the  user's  program  as: 

DEFINE  CALL  &NAME  (&PARAMETERS)  = 
&IF  &SIGN  (&MFIELD(&NAME)) 

&THEN  EXTERNAL  &NAME;  &FI 
&IF  &EMPTY  (&PARAMETERS)  &TKEN  &ELSE 
BEGIN  BLOCK 

BEGIN  USE  (63) 
LIST:   DATA  & PARAMETERS 
END; 
CLC(2); 
SLIT (2)  LIST; 
END;  &FI 
CLC(3); 

SLIT (3)  &NAME; 
EXCHL(3)  &ICR  ##; 

2.2  Core  storage  used 

This  routine  uses  500  words  of  PE  memory  for  storing  instructions. 
One  row  is  used  for  storing  PE  numbers  and  three  additional  rows  are  used 
for  temporary  storage.   ADBO  ~  32  are  also  used. 

2.3  Constant 

EPS  (e),  the  constant  which  is  used  to  test  the  convergence  (k.l), 
is  taken  as  10"   in  this  program.   If  a  user  wants  to  change  the  value  of  this 
constant,  he  may  insert  EPS:   DATA  (desirable  value);  in  place  of 

EPS:   DATA  @  -  10;  . 


3.   EXAMPLE 


matrix: 


A  test  of  this  program,  was  made  on  the  B5500  simulator  for  the 


5-0 

-2.0 

-5.0 

-1.0 

1.0 

0.0 

-3-0 

2.0 

0.0 

2.0 

2.0 

-3-0 

0.0 

0.0 

1.0 

-2.0 

with  €  =  10   .   The  comparison  of  the  eigenvalues  obtained  by  this  program  to 
the  exact  values  is  given  in  Table  1.   The  selection  of  this  small  matrix  and 
a  relatively  large  e  was  made  because  of  the  speed  of  the  SSK  simulator  on 
the  B5500;  the  execution  speed  ratio  of  the  simulator  to  the  ILLIAC  IV  is 
approximately  1:10  . 

Table  1 


Eignevalues  obtained 

Exact  eigenvalues 

on  B5500  simulator 

3-999999997867 

k.O 

I.OOOOOOOOIO66  +  1. 99999999992 8i 

1.0  +  2.0i 

I.OOOOOOOOIO66  -  1. 99999999992 8i 

1.0  -  2.0i 

-1.000000000000 

-1.0 

h.      QR- ALGORITHM 

^.1  Brief  outline  of  the  QR-algorithm 

The  QR-transformation  consists  of  the  decomposition  of  the  matrix 
A  into  the  product  of  a  unitary  matrix  Q  and  an  upper  triangular  matrix  R  , 
and  forming  A      by  post-multiplying  R  by  &.   Thus 

\  +  i  -  W       where      \  ■  \\'  (1) 

therefore 

where  A  =  A  is  the  original  matrix.   It  can  be  shown  in  general  that  A 

(k)       (k) 

tends  to  a  form  in  which  a.  '        .    a.  '        .        n  =  0  for  i  =  0,  1, . . . . ,  N  -  3 

as  k  increases.   All  eigenvalues  are  therefore  either  isolated  on  the  diagonal 
or  they  are  eigenvalues  of  a  2  x  2  diagonal  submatrix.   The  amount  of  calcula- 
tions involved  in  a  QR  step  is  greatly  reduced  if  the  matrix  A  is  in  the 
Hessenberg  (or  almost  triangular)  form.   Since  there  are  several  stable 
methods  available  to  reduce  a  general  matrix  to  this  form  (ASK  program.  HSBG 
is  written  for  this  purpose),  the  QR-algorithm.  is  used  after  such  reduction. 

In  order  to  achieve  rapid  convergence,  it  is  essential  that  the 
origin  shifts  be  applied  and  that  each  shift  be  close  to  an  eigenvalue 
of  the  matrix.   The  QR-algorithm  with  shift  of  an  origin  s  is  expressed 
as: 

\  + 1  -  \  \  +  \z      where      \  -  V  =  W  <2> 

or  in  other  words 

\  +  i  ■  \  \V 

However,  even  when  A  is  real,  some  of  the  eigenvalues  may  be 

complex.   If  the  transformation  (2)  is  carried  out  with  a  complex  value  of 

s  ,  A     is  in  general  a  complex  matrix.   This  deficiency  can  be  overcome 

by  performing  two  steps  of  (2)  with  shifts  of  s  and  s      respectively. 

Since  s.  and  s,    .  are  both  real  or  complex  conjugate  in  this  transformation, 
k      k  +  1  7 

A     should  be  always  real.   This  transformation  is  described  as 


and 


<\\  ♦  1>  <Rk  +  1  V  =  <\  "  S?>  <\  -  sk  + 11)  (3) 


One  method  of  calculating  A^     by  (3)  would  be  to  form  the  real 
matrix  r  =  (A  -  si)  (A,  -  s     i),  computing  its  unitary- triangular 
decompositon  to  obtain  Q,  Q      and  transform  A,  by  means  of  this,  thus 
giving  A^   p.   This  process  requires  a  prohibitive  amount  of  work,  but  it  is 
shown  [l]  that  when  the  matrix  is  in  the  Hessenberg  form,  it  is  unnecessary 
to  compute  more  than  the  first  column  of  r,  and  that  this  immediately  gives 
the  transformation  to  be  applied  to  A^. 

k.2     Practical  computation 


If  (3)  is  rewritten  as 

\  +  2  =   w*Akw     and     W*r  =  A 

where  W  =  Q,  Q.    .  and  A  is  the  triangular  matrix  R,     R  ,  W  is  a  unitary 
matrix  which  reduced  r  to  the  triangular  A,  and  W  is  composed  of  N  unitary 

Fi   1  °  " 

factors  of  the  form  M.  =  so  that  W  =  M,  Mp  M.  .   From  the 

L        i-" 

form  of  each  M.  we  see  that  the  first  column  of  W  is  equal  to  the  first 

column  of  IVL ,  and  this  is  any  unitary  matrix,  the  transpose  of  which 
eliminates  the  elements  of  the  first  column  of  V   below  diagonal. 

Since  we  wish  to  transform  A^  to  A      by  W,  we  first  operate  on 
A  with  M  .   This  will  change  the  first  three  rows  and  columns  of  A^  since 
the  first  column  of  r  contains  only  three  non-zero  elements.   It  follows; 


V 


00 


"10 


01 


11 


21 


l02  * 


'12 


22  ' 


.  .  a 


0,N  -  1 

1,N  -  1 

l2,N  -  1 


'N-  1,  N  -1 


-» 


00 


01 


02 


aio 

!ii 

a!2 

a20 

a21 

tk 

a30 

a31 

a32 

"03 


°13 


£3 
"33 
l43 


0,N 
ll,N 
l2,N 


3,N  - 


N  -  1,  N  -  1 


where  the  elements  changed  by  the  row  and  column  operations  are  underlined 
and  primed  respectively.   The  resulting  matrix  is  no  longer  in  the  Hessenberg 


form  and  since 


\ 


+  2 


is  in  the  Hessenberg  form,  we  can  say  that  the  matrices 


M2.  M3 


.  .  NL  reduce  M.  A^  M  to  the  Hessenberg  form.   In  the  practical 
computation,  therefore,  it  is  necessary  to  compute  only  the  first  column  of  r. 

After  each  iteration  (we  call  the  calculation  of  A,    _  from  A 

v  k  +  2      ± 

an  iteration),  all  subdiagonal  elements  of  A,    p  are  examined  to  see  if  any 
of  them  are  "negligibly  small."  If  so,  the  eigenproblem.  for  the  current 
matrix  splits  into  that  for  two  or  more  Hessenberg  matrices  of  smaller  sizes, 
and  the  iterations  continue  with  the  submatrix  in  the  bottom  right-hand 
corner.   It  may  happen  that  while  no  individual  subdiagonal  element  is 
sufficiently  small  to  be  regarded  as  negligible,  the  product  of  two  consecu- 
tive elements  may  be  small  enough  to  permit  us  to  work  with  a  submatrix. 
Therefore  the  examination  of  the  matrix  A     is  performed  to  see  if  any 
two  consecutive  subdiagonal  elements  are  small. 


5-   PROGRAM  DESCRIPTION 

5.1  Search  for  negligible  subdiagonal  elements 

We  assume  that  the  size  of  the  matrix  under  consideration  is 
(n  +  l)  x  (n  +  l)  where  n  takes  integer  values  between  1  and  N  -  1.   If  the 
last  negligible  subdiagonal  element  is  in  position  (£,    £  -   l),  it  is  required 
only  to  work  on  the  submatrix  in  the  rows  and  columns  £   to  n.   If  none  of  the 
subdiagonal  elements  are  negligible,  I   is  taken  to  be  0.   The  following 
criterion  is  used 

|aJ,  J-il  <  e^ai  -  l,  i-  J  +  lag,  jl); 
This  criterian  examines  whether  a,£     &   _   is  negligible  compared  to  the  local 
diagonal  elements. 

On  each  PE  £,    the  following  computations  are  simultaneously 
performed: 

f(i)  S-|aif  £_±\    -     e  (\,£_   ^   £_   1|  +  |a^  £\ 

for  i  =  1,  2,    .    .  .  ,  n. 

If  f(i)  is  negative,  1  is  placed  in  ith  bit  of  the  ACAR.   Then 
searching  is  made  for  the  lowest  bit  of  the  ACAR  which  contains  1.   For 
example,  in  the  following  case: 

0   1   2   3  ...  19   20 63 

ACAR 


0   0   1   0...0    1   0   0...0   0 


I   is  set  to  be  20. 


Then  test  is  made  if  £  =   n  or  Z  =   n  -  1.   If  £  =   n,  one  eigenvalue 
is  found  in  the  place  (n,  n)  and  the  matrix  is  deflated  by  1,  and  n  is 
decreased  by  one.   In  the  case  that  £  =   n  -  1,  two  eigenvalues  are  found  as 
the  eigenvalues  of  the  bottom-right  hand  corner  2x2  submatrix.   Then  two 
columns  and  rows  are  deflated  and  n  is  decreased  by  2. 


r> 


5.2  Shifts  of  origin 

The  shifts  of  origin  at  each  stage  are  taken  to  be  the  two  roots, 

s.  and  s,    n ,  of  the  2x2  matrix  in  the  bottom  right-hand  corner  of  the 
kk+1  ° 


current  iL.   This  gives 


s.  +  s.    t  =  a 

k    k  +  1    n 


1,  n  -  1 


+  a 


n,  n 


and 


s,  s.    _  =  a 

kk+1    n 


1  _,  a   -  a    n    a 

1,  n  -  1  n,n    n  -  1,  n  n,  n-1 


(h) 


In  some  rare  cases,  the  process  fails  to  converge  with  these  shifts 
of  origin.   An  example  of  such  failure  is  provided  by  matrices  of  the  type: 


0 

0 

0 

0 

1  " 

1 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

1 

0 
0 

Here  the  shifts  of  origin,  given  in  (k) ,  are  both  zero,  and  since  the  matrix 
is  orthogonal,  it  is  invariant  with  respect  to  the  QR  transformation  without 
shifts.   However,  if  one  iteration  is  performed  with  any  shifts  of  origin 
which  are  loosely  related  to  the  norm,  of  the  matrix,  the  convergence  is  very 
rapid.   Therefore  in  the  case  where  ten  iterations  do  not  produce  an  eigen- 
value, the  usual  shifts  s.  and  s,    _  are  replaced  by  shifts  defined  by 
'  k      k  +  1 


s.  +  s.    .,  =  1.5  (  a       n   +  a    _ 
k    k+1       v  '  n,  n-11    'n-1,  n 

2 


2 


)  • 


s,  s. 


=  (U. 


I) 


(5) 


k  k  +  1   Vl  n,  n  -  l1    '  n  -  1,  n  -  21 
This  strategy  is  used  again  after  20  unsuccessful  iterations.   If  30 
unsuccessful  iterations  are  needed  then  a  failure  indication  is  given. 

In  this  program,  ITS  is  the  name  of  an  ADB  where  the  iteration 

count  is  stored.   When  ITS  +-  10,  20,  we  form  S  =  s.  +  s.    ,  and  Y  =  s.  s,    . 

'  kk+1         kk+1 

on  a  PE  n  according  to  (k)   and  store  in  ADB's.   When  ITS  =  10  or  20,  scheme 
(5)  is  used  for  computing  S  and  Y.   If  ITS  =  30,  it  is  assumed  that  this 
algorithm  fails  to  produce  eigenvalues.   As  a  result,  only  eigenvalues 
computed  prior  to  this  point  are  given  in  WR  and  WI. 


10 

5-3  Search  for  two  consecutive  small  sub diagonal  elements 

After  determining  &,    (5-l),  the  submatrix  in  the  rows  1  to  n  are 
examined  to  see  if  any  two  consecutive  subdiagonal  elements  are  small  enough 
to  work  with  an  even  smaller  submatrix.   To  test  if  we  are  to  start  at  the 
row  m,  we  compute  the  elements  p  ,    q  and  r  such  that 

p=a   -a   {s.    +  s,    _ ;  +  s.  s,    ,  +  a       a 
m    mm.    mm   k    k  +  r  k  k  +  1    m,  m  +  1  m  +  1,  m 

q.  =  a        (a   +  a    _      ,  -  s,  -  s,    n  )  (6) 

m    m  +  1,  m  v  mm    m  +  1,  m  +  1    k    k  +  ly  K    ' 

m"      m  +  2 ,   m  +   1     m  +   1,   m.   " 
The  criterion  applied  is 

i  I     C  f  <3 I    +    lr    I  ) 

-  1  tti  '    m ' 


m,    m 


(7) 


<    €    |p   |    (|a  |    +    la        !   +    la 

—  Irm'  m  +   1,   m  +   l1         '    m,m'         '    m  -   1,   m  -l|  ) 

where  we  test  whether  or  not  the  elements  which  appear  in  the  positions 

(m.  +  1,  id),  (m  +  2,   m  +  l)  are  negligible  compared  with  the  three  local 

diagonal  elements  a    _      _ ,  a    and  a 

m  +  1,  m  +  1'  m,m      m  -  1,   m  -  1 

Here  we  take  m  to  be  the  largest  integer  (>  i)  for  which  condition  (6)  is 

satisfied. 

For  this  computation,  the  mode  bits  are  turned  on  for  PE  I   through 
PE  n.   The  pm,  q^  and  r  for  i  =  1,  I  -  1,  .  .  .  ,  n  are  computed  according 
to  (6)  on  all  PE's  whose  mode  bits  are  turned  on,  and  comparison  is  made  to  see 
whether  (7)  is  satisfied.   If  (7)  is  satisfied  on  PE  m,  1  is  placed  in  the 
mth  bit  of  the  ACAR.   Then  the  search  is  made  for  the  lowest  bit  of  the 
ACAR  which  contains  1,  and  m  is  set  equal  to  this  bit  number.   If  no  1  is 
found  in  the  ACAR,  m.  is  taken  to  be  i. 


11 

5.+     Double  QR-transformation 

A  is  computed  by  applying  the  QR-double  transformation  to  iL 

in  such  a  way  that  A.         0=N     .    .    .    .   N    A,    N     .    ...    IT 

k  +  2         n  m     k     m  n 

* 

*  UiUi 

where  N.    =  I ?r  and  U.    =   (p.    +  t .  ,    q.  ,   r .  ,    0   .    .    .    .0). 

l  0„  2  l         XJri  —     i;    ^i'      i' 

1 


Here 


and 


p.  =  a.  .    -   a.  .    (s.    +  s.         . )  +  s,    s,         _   +  a.    .        .    a. 
i         n         ii       k         k  +   1  k     k  +   1         1,1  +1     i  +  1,: 

q.  =  a.        _    .    (a.  .    +  a.        _    .        .    -   s,    -  s,         n) 
*i         i  +  l,i   v  ii         i+l,i+l         k         k  +  r 

r.  =  a.        _    .         _    a.         _    .  for   i   =  m 

l  i  +  2,i  +   l     i  +   l,i 


p.    =  a.    .             +  t. ,  q.    =  a.    ,           . 

i          i,i-l—i'  l          i  +   1,    l   -   1 

and  r.  =  a.        _    .        _                     for   i  *  m. 

l  l  +  2,i   -    1                                  T 


t.    and  2K.    are  defined  as 
i  i 


/2  2  2 

t.    =  +    ./p.    +   q.    +   r. 
x        —    v      J.  1  1 

2  2  — 

2KT  =  t.    +  p.t. . 
i  ill 

Row  modification: 

*   (i) 
For  i  =  m,  m  +  1,  .  .  .  ,  n,  the  elements  of  N.  A  .  = 
'  '  '     '  i   k 

N*  (N.    ,  •  .  .  N*  A,  N  .  .  .  N.    ,)  are  different  from  those  of  A\X'  in 
ii-l       m  k  m       l-l  k 

only  three  rows,  i.e.,  ith,  (i  +  l)th  and  (i  +  2)th  rows.   These  new  elements 
are  computed  in  the  following  way  with  the  elements  of  A    denoted  by  a,  . : 

(i, j )  -  element  =  a. .  -  [(p.  +t.)a..  +q.a.  .  .  +r.a.  _  .]— 


(i+l.j  )  -  element  =a.  n  .  -  [(p.  +t.)a..+q.a.n  .  +  r.a.  _  .]  ^ 

V  '°  1+1,0  1    ~        1         1J  1     1+1,J  1     1+2,  J         pir     2 


;i+2,j)    -   element   =  a.+^  .    -    [(p.    +  t.Ja..   +  q.a.+^  .   +  r.a.+^  .}  _fi 

2K. 

l 

for   j   =   i,i  +   1,    ...,n 
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2 
In  the  actual  computation,  t.  and  2K.   are  first  computed  and  then 

p  ,  q.  and  r.  are  found  and  stored  in  ADB's.   The  mode  bits  of  PE's  which 
*V   1      l 

contain  a.  .,a.  .  .,,  .  .  .  .  ,a.   are  turned  on  and 
1,1   i,i+l  i,n 

c.  :  (p.  +  t. )a. .  +  q.a.  _  .  +  r.a.  _  .  are  computed  on  these  PE's.   The 
l  -  v  i  -  i  13    i  i+l,  3    i  i+2,J 

computation  of  new  elements 

c. 

(i.j)  -  element  =  a.  .  -  t— 
v  '  ij   t. 

e.g. 

-i  -*"i 

(i.j+l)  -  element  =  a.  .  ,.  -  ?r 

1>J  +  1   2K 
i 

c.r. 
(i,j+2)  -  element  =  a.±    .  2 ^ 

2K . 

l 

are  then  performed. 

Column  Modification: 

/  *    (Ox  *   (i) 

Similarly  (N.  A  v   )N.  is  computed  from  N.A  v  '   for  l  =  m,  m  +  1, 

*   (k) 
.  .  .  ,  n  in  the  following  way  where  the  element  of  matrix  N.A.     are 

denoted  as  a .  ,  : 
3,h 

(d,i)   -   element   =  a..    -   [  (p.   +  t.)a..    +  Vj,i+1  +  riaj,i+2]  \ 

qi 

(j.i+l)   -   element   =  a.    .    ,    -    [(p.    +  t.)a..    +  q.a.  +  r.a.  ]  5 

Vd'        '  3,i+l  v*i  -     i'    ji       *i  3,i+l         i   3,i+2     2K  2 

i 


(3,1+2)  -  element  =  a     -  [  (p.  +  t.)a..  +  q^^  +  r^^]  -ig 


for  j  =  i,    .  .  .  •  ,  min  [i+3,n]  . 

As  in  the  row  modification,  c!  =  (p.  +  t. )a. .  +  q.a  .  .  +  r.a      are 

l     i  —  l  31    1  3,1+1    1  3,1+^ 

computed  on  the  PE's  which  contain  a.,  for  3  from  i  through  min  [i  +  3,  n], 

then  the  computations  of  new  (3*,i)  -  element,  (3,  i  +  l)  -  element  and 

(j,  i  +  2)  -  element  are  performed  on  the  corresponding  PE's. 

5 • 5  Computation  of  eigenvalues 

The  eigenvalues  are  calculated  as  the  last  step  of  program  after 

a.  ,  .  or  a.  .  ,  a.  ,  .  become  negligibly  small  for  all  0<i<2i-2.  At  each 
1+1,1     i,i-l  1+1,1         to 
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time  when  the  matrix  is  deflated  by  2,  1  is  placed  in  position  n  of  ADB  named 
SOLV.  If  the  transformation  is  carried  out  successfully,  SOLV  looks  like  the 
following, 

0123^56789  10  62  63 


SOLV 


001000010010 0      1 


Here,    the   eigenvalues  we  are   looking  for  are  a^_,    a__,    a,  ,  ,    a,.,., 
t  &  &  qq>     22>     4V     55' 


a88>  all,ll' 


,  a,-,  £,  and  eigenvalues  of  the  following  2x2  matrices 
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"22 


%6 


77 


'99 


10, 10 


and 


a 


'62,62 


l63,63 


These  eigenvalues  can  be  calculated  simultaneously  on  the  corresponding  PE's. 


5.6  The  number  of  iteration  required 

After  each  deflation,  the  number  of  iterations  required  to  find  each 
eigenvalue  (or  two  eigenvalues)  appears  in  ADB  named  ITS.   This  is  placed  as 
the  nth  element  of  the  row  vector  specified  by  the  user  as  the  fifth  para- 
meter in  the  calling  sequence.   If  this  parameter  is  specified  as  0,  the 
number  of  iterations  is  not  stored  anywhere. 
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n  «-  size  of  matrix  -  1 
clear  rows  IT  and  SOLV 


Yes 


ITS  *-  0 


look  for  single 
small  subdiagonal 
element   and  find  I 


Yes 


Yes 


Yes 


FIN, 


ONEW] 


TWOWJ 


A  I] 


Yes 


No 


computer  origin  shifts 

V  sk+i  and 

s  =  sk +  sk+l 
Y  =  sk  sk+l 


compute  origin  shifts 

S  =  L'5  (l  an,n-ll  +  I  an-l,n-2') 

Y;:  Cla^   x|  +|an.^n.2D2 


ITS  <-  ITS  +  1 
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i  «-  i  +  1 


no 


1 


look  for  two  consecutive 
small  sub diagonal 
elements  and  find  m 


i  <-  m 


yes 


no 


p,  q,  r  were  already 
computed  in 

FINDM 


compute  p,  q  and  r 


compute  t  and  2K 


_v- 


row  modification 
compute  new  element 
a .  . ,  a.  n    .,  a.  _ 

io'  i+i, a'  1+2, j 

for  j  =  1,    .  .    .  ,  i 


^ 


column  modification 
compute  new  element 

a.  ..  a.  .   «  a.  .  ~ 

D,i   D,!4"!   a,  1+2 

for  j  =  i,  .  •  •  ,  n 


IT 


(oimT) 

-^  fifth  argument 

calling 


no 


store  ITS  into  row  IT 
as  an  nth  element 


±k- 


n  <-  n  -  1  <; 


(  NEXTW  J 


yes 


-* 


yes 


no 


store  ITS  into  row  IT 
as  an  nth  element 


store  1  into  nth  bit 

of  SOLV 


ik_ 


n  <-n  -  2 


C  NEXTWJ 
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FIN 


mode  bit  <-  1  for  all  PE's 
■which  satisfy 

0  <  PEN  <  N  -    1 


find  eigenvalues  on 
all  enabled  PE's 


WR  <~  real  part  of 
eigenvalues 

WI  «-  imaginary  part 
of  eigenvalues 


FAIL 


mode  bit  «-  1  for  all  PE's 
which  satisfy 

n  <  PEN  <  N  -  1 


end  of  this 
subroutine 
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